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We study the phase diagram of calcium fluoride (CaF2) under pressure using classical molecular 
dynamic simulations performed with a simple pairwise interatomic potential of the Born-Mayer- 
Huggings form. Our results obtained under conditions < P < 20 GPa and < T < 4000 K 
reveal a rich variety of multi-phase boundaries involving difi^erent crystal, superionic and liquid 
phases, for all which we provide an accurate parametrization. Interestingly, we predict the existence 
of three special triple points (i.e. solid-solid-superionic, solid-superionic-superionic and superionic- 
superionic-liquid coexisting states) within a narrow and experimentally accessible thermodynamic 
range of 6 < P < 8 GPa and 1500 <T < 2750 K. Also, we examine the role of short-ranged repulsive 
(SR) and long-ranged attractive (LA) atomic interactions in the prediction of melting lines with the 
finding that SR Ca-F and LA F-F contributions are most decisive. 

PACS numbers: 66.30.H-, 81.30.Dz, 62.50.-p, 45.10.-b 



I. INTRODUCTION 



Calcium fluoride (CaF2) is representative of the 
fluoride-structured halides, an important family of ionic 
materials with numerous applications in high-pressure 
science and technology. Examples of its assorted range 
of qualities include, excellent optical transmission proper- 
ties over a wide wavelength range, a large electronic band 
gap, and very high elastic compressibility [iHll- CaF2 is 
also well known for being a fast-ion conductor, a material 
in which the lighter ions (i.e. F~) acquire a significant 
mobility comparable to ionic melts at temperatures well 
below its fusion point [^-lll- This unusually large dif- 
fusion of fiuor anions through the almost rigid matrix of 
calcium cations (i.e. Ca+^), an effect also termed as supe- 
rionicity, is originated by the formation of Frenkel pairs 
and migration of interstitial ions (rl-fioj. Superionicity, 
which may be also found in other oxide- , hydride- and 
iodide-based complexes like Y2O3, LiBH4 and Agl [TT| - 

, find an use in solid state applications as supercapac- 
itors, batteries and fuel cells. 

The intrinsic variety of condensed matter phases in 
CaF2 (i.e. solid, superionic and liquid) a priori sug- 
gests a rich and very intriguing P ~ T phase diagram. 
Nevertheless, most of the experimental and theoretical 
investigations performed so far have focused on very nar- 
row low-pressure and low-temperature thermodynamic 
ranges. The T = polymorphism of CaF2 encompasses 
two main crystal structures (see Fig.[T]), the low-pressure 
face-centered cubic fiuoride phase (space group Fmim) 
and the high-pressure orthorhombic PbCl2-type (cotun- 
nite) phase (space group Pnma). Although CaF2 fast- ion 
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conduction has been intensively investigated at ambient 
pressures (6l-[To|. it surprisingly remains yet unexplored 
in the cubic and orthorhombic phases under compres- 
sion. Gathering information on the pressure dependence 
of CaF2 fast-ion conduction and melting however turns 
out to be very desirable both from a fundamental and 
technological point of view. Assuming fast-ion conduc- 
tion behavior in both low- and high-pressure phases, for 
instance, opens the possibility for the existence of spe- 
cial triple or even quadruple points (i.e. thermodynamic 
states in which several superionic and solid phases might 
coexist in thermodynamic equilibrium) at elevated P — T 
conditions. Also, our present knowledge on the fusion 
properties of CaF2 and of ionic compounds in general, 
is unfortunately rather scarce Likewise, the ex- 

pected connections between atomic sub-lattice melting 
and homogeneous melting, a subject that we have re- 
cently investigated in Ar(H2)2 [3, are not yet clearly un- 
derstood. On the practical side, pressure-induced trends 
unravelled in CaF2 could be generalized to other fiuoride- 
structured materials composed of heavy and light ion 
species, that have been proposed or are used in relevant 
technological applications (e.g. transition metal hydrides 
in hydrogen storage cells and uranium dioxide for gener- 
ation of nuclear fuel). 

In this article, we report the phase diagram of CaF2 
under pressure (e.g. < P < 20 GPa) as obtained 
from classical molecular dynamics (MD) simulations per- 
formed with a simple but reliable rigid-ion interatomic 
potential of the Born-Mayer-Huggings (BMH) form. Our 
results rely on extensive one- and two-phase coexis- 
tence simulations and reveal a rich variety of p revi- 
ously unreported P — T phase boundaries [l5l - [l7l | (i.e. 
fiuoride-PbCl2, fluoride-liquid, PbCl2-liquid, fiuoride- 
superionic fiuoride, PbCl2-superionic PbCl2, superi- 
onic fluoride-PbCl2 and superionic fluoride-superionic 
PbCl2), for all which we provide here an accurate 
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Fluoride 
(Fm3m) 




PbCI^-type 
(Pnma) 





A (eV) 


P(A) 


C (eV A*') 


Ca-F 1717.441 


0.287 


0.102 


F - F 2058.994 


0.252 


16.703 



Ca 



TABLE I: Interatomic pair potential parameters for 
CaFz [13. 



next section, we describe the computational methods em- 
ployed and the low-temperature performance of the se- 
lected Born-Mayer-Huggings potential. In Sec. IIIIl we 
present our results for the phase diagram of CaF2 un- 
der pressure and discuss them. Then, a section follows 
in which we analyze the role of repulsive and dispersion 
interactions on the determination of melting points by 
theoretical means. Finally, we summarize the main con- 
clusions in Sec. |Vl 



FIG. 1: Cubic fluoride and orthorhombic PbCl2-type (con- 
tunnite) structures of CaF2. In the fluoride structure Ca^^ 
ions are in a cubic close-packed arrangement and F~ ions 
occupy all the tetrahedral interstitial sites. In the contun- 
nite structure the array of calcium cations is hexagonal close- 
packed (hep) and half the fluor anions is placed off-center 
in the ideally octahedral hep voids with fivefold coordination 
(the other half of fluoride ions exhibiting tetrahedral coor- 
dination). In the fluoride structure Ca atoms are eightfold 
coordinated whereas in the contunnite structure the coordi- 
nation number increases to nine. Calcium and fiuor ions are 
represented with blue and grey spheres respectively, and the 
unit cell with black solid lines. 



parametrization. For example, we find that both fluoride- 
superionic fluoride and PbCb-superionic PbCl2 phase 
boundaries, Ts{P), are linearly dependent on pressure 
with a small and positive dTg/dP slope of 34.2 and 
50.2 K/GPa, respectively. Interestingly, we predict the 
existence of three special triple points involving coex- 
istence of fluoride-superionic fluoride-PbCl2, superionic 
fluoride-PbCl2-superionic PbCl2 and superionic fluoride- 
superionic PbCl2-liquid phases, within a narrow and ex- 
perimentally accessible region of 6 < P < 8 GPa and 
1500 <T< 2750 K. In addition to these findings, we use 
a free-energy perturbative approach to evaluate the shift 
in melting temperature caused by mild variations of the 
employed potential parameters. By doing this, we quan- 
tify the role of repulsive and dispersion interactions on 
our results and qualitatively gain access to the melting 
features of BMH potentials used in other works. 

The organization of this article is as follows. In the 



II. SIMULATION DETAILS AND THE 
INTERATOMIC PAIR POTENTIAL 

Calculations were done with LAMMPS [l^, a paral- 
lel classical molecular dynamics (MD) code comprising a 
large variety of potentials and different schemes for sim- 
ulation of solid-state and soft materials. Our MD cal- 
culations are of two main types: one-phase (i.e. pure 
liquid, solid and superionic phases) and two-phase co- 
existence (i.e. liquid and superionic phases coexisting 
in thermodynamic equilibrium) simulations. One-phase 
simulations are performed in the canonical {N, V, T) en- 
semble and two-phase coexistence simulations in the mi- 
crocanonical {N, V, E) ensemble (specific details of these 
simulations are provided in Sec. IIIip . In [N^V^T] simu- 
lations the temperature is kept fluctuating around a con- 
strained value by using Nose-Hoover thermostats. Large 
simulation boxes containing 6, 144 and 12, 288 atoms are 
used in our one-phase and two-phase coexistence sim- 
ulations, respectively. Periodic boundary conditions are 
applied along the three Cartesian directions in all the cal- 
culations. Newton's equations of motion are integrated 
using the customary Ver let's algorithm and a time-step 
length of lO^'^ ps. A particle-particle particle-mesh k- 
space solver is used to compute long-range van der Waals 
and Coulomb interactions beyond a cut-off distance of 
12 A at each time step. 



(1) 



The interatomic potential adopted for this study is 
U{r) = VcaF(?') + ^ff('') where terms Vij are of the 
Born-Mayer-Huggings (BMH) form (see Eq. [T|). Each 
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pairwise term is composed of three different contribu- 
tions; the first is of exponential type and accounts for the 
short-ranged atomic repulsion deriving from the overlap- 
ping between different electron clouds; the second term 
is proportional to r~^, with r being the radial distance 
between a given couple of ions, and represents the long- 
ranged atomic attraction due to dispersive van der Waals 
forces; the third term is the usual Coulomb interaction 
between puntual atomic charges, which in our case are 
taken to be = +2e and Zp = — le. In Table I, we 
enclose the value of the BMH parameters used through- 
out this work and which coincide with those primarily 
deduced by Dick and Overhauser [2^. It must be noted 
that the original Dick- Overhauser potential includes elec- 
tronic polarization effects via a shell model however for 
present purposes these can be safely neglected due to the 
marked ionic nature of CaF2. Indeed, Lindan and Gillan 
explicitly showed for a similar BMH model that inclusion 
of electronic polarizability had a remarkable small effect 
on the estimation of static and dynamic CaF2 quanti- 
ties 0,11. 

In order to assess the reliability of the adopted BMH 
potential, we performed a series of static ground-state 
calculations and compared them with available low-T ex- 
perimental data. In Fig. [21 we show the zero-temperature 
equation of state of CaF2 obtained for its cubic and or- 
thorhombic phases (solid lines). In each phase, we com- 
puted the energy per formula unit for a set of 20 volume 
points spanning over the range 9.0 < < 14.5 A^. Sub- 
sequently, we fitted the results to a third order Birch- 
Murnaghan equation of the form 



^pcrf(l^) = Eo + -VaKo ■ 



2(1 + X) y 



(2/3) 



+ 



1 



x + 



(2) 



where Eq and Kq = ^Vq^^ are the energy and bulk 
modulus at equilibrium volume Vo, x 



and 



Kq = [dK/dP], with derivatives taken at zero pressure. 
(Atomic forces and cell shape relaxations were performed 
for the orthorhombic PbCl2-type phase at each volume.) 
The value of the resulting Eq, Vq, Kq and Kq parameters 
are -9.05 (-8.97) eV, 13.56 (12.52) A^, 108.1 (79.5) GPa, 
and 1.51 (5.81) for the cubic (orthorhombic) structure. 
The static equation of state then is obtained as the mi- 
nus derivative of Eq. ([2]). By comparing the enthalpy of 
the different phases we find that the zero-temperature 
cubic — > orthorhombic phase transition occurs at a pres- 
sure of Pt = 10.95 GPa, in very good agreement with 
recent ambient experimental data obtained by Kavner 
and others [U, [2^ . The corresponding change of volume 
is -8.35 % with V = 12.33 (11.30) A^ in the cubic (or- 
thorhombic) phase. 
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FIG. 2: Zero-temperature equation of state of CaF2 un- 
der pressure obtained with the present BMH potential (solid 
lines). Experimental data can be found in Refs [lH] and [2^ . 



In Table II, we report the bulk modulus and struc- 
tural parameters of the two polymorphisms of CaF2 ob- 
tained at equilibrium and Pt- As one can see, the ac- 
cordance with experiments in this case is also notable. It 
must be noted that from the Birch-Murnaghan fit quoted 
above we obtain a zero-pressure bulk modulus that is 
~ 20 GPa larger than the experimental value reported 
by Dorfman et al. [13]. However, in the obtaining of that 
experimental datum the value of Kq parameter in the 
corresponding equation of state fit was set to 4.7. Pro- 
ceeding in the same way as Dorfman and collaborators, 
we obtain Kq = 82 GPa which is in very close agree- 
ment with their measurements. Finally, we computed 
the value of the three independent elastic constants of 
the Fm3m phase at equilibrium, C^j's. For this, we dis- 
torted the shape of the cubic unit cell according to the 
strain matrices shown in Table HI and fitted the result- 
ing variation of the energy to the also reported parabolic 
curves We obtain Cn = 168.1, C12 = 46.9, 

and C44 = 37.7 GPa which agree notably with the ex- 



44.4, and 



perimental values C^i^* = 165.4, C- 
C|7* = 34.2 GPa [H, and previous theoretical estima- 



expt 
12 



tions as well [17 1 



The main conclusion emerging from these calculations 
is that the adopted BMH potential provides a very reli- 
able account of the low-T region of the phase diagram of 
CaF2. It can then be assumed that medium and high- 
temperature descriptions obtained with this same simple 
model interaction will be also physically meaningful. In 
fact, as we will show in the next section, superionic and 
melting temperatures predicted at ambient pressures are 
fully consistent with experimental observations. 
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Theory 


Experiment 


Theory 


Experiment 


FmSm (P = 


--0) 


ao = 5.460 


ao = 5.463" 


108 (7^0 = 1-5) 
82 {K'o = 4.7) 


85' {k'o = 4.7) 


FmSm (P = 


Pt) 


ao = 5.289 


ao = 5.313'' 


122 


152(20)== 


Pnma (P — 


Pt) 


ao = 5.721 

bo = 3.463 

Co = 6.846 
Ca (0.2472,0.25,0.1153) 
Fl (0.8510,0.25,0.0732) 
F2 (0.4768, 0.25, 0.8300) 


ao = 5.700'' 

bo = 3.450'' 

Co = 6.800'' 
Ca (0.2530, 0.25, 0.1094)^* 
Fl (0.8595,0.25,0.0731)'* 
F2 (0.4780, 0.25, 0.8344)=* 


138 


162(30)== 



TABLE II: Structural parameters and bulk modulus K — —V{dP/dV)T of CaF2 obtained at zero-temperature; theory values 
are obtained with the present BMH potential and experimental uncertainties are indicated within parentheses. Pt = 10.95 GPa 
is the pressure at which the cubic — > orthorhombic phase transformation is found to occur. Experimental data can be found in 
Refs. [ill", [23^ and [H'^. 
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u{e) = I (Cll +2Ci2)e2 



u(e) = 3 (Cll - C12) e^ 



M(e) = 2C44e^ 



TABLE III: Strain matrices for the three independent elastic constants Cn, 
corresponding strain energy relations per unit volume at P = 0. 



C12 and C44 of the cubic Fm3m structure and 



III. RESULTS AND DISCUSSION 

In the next sections we will describe in detail the two- 
phase boundaries and special triple points appearing on 
the phase diagram of CaF2 shown in Fig. [31 however let 
us first explain the general strategy that we followed to 
obtain it. 

Initially, we performed extensive one-phase and two- 
phase molecular dynamics simulations to find out the 
solid-superionic and superionic-liquid phase boundaries 
of the fluoride and PbCl2-type structures in the whole 
range of pressures considered (i.e < P < 20 GPa). 



Next, we considered the thermodynamic state (Pic,T'ic) 
at which the melting curves of the two superionic states 
cross each other. The Gibbs free energy of the supe- 
rionic fluoride and superionic orthorhombic phases are 
equal at that thermodynamic state thus (Pic, Tic) must 
belong also to the boundary separating the superionic 
cubic and superionic orthorhombic regions. In other 
words, (-Picnic) is a special triple point {special be- 
cause it comprises coexistence of superionic and liquid 
phases, in contrast to ordinary solid-solid-liquid and/or 
solid- liquid- vapor triple points). One-phase MD simula- 
tions were then carried out to determine the volume and 
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FIG. 3: The phase diagram of CaF2 under pressure obtained from molecular dynamics simulations. Solid lines represent 
function fits (see text) and filled dots thermodynamic states at which molecular dynamic simulations were explicitly carried 
out. Transition temperature errors are represented with solid bars. 



enthalpy (i.e. H = E + PV) of the two superionie phases 
at (Pic, Tic). Using the customary Clausius-Clapeyron 
relation 



dT 
dP 



T 



AV 
AH 



(3) 



we computed the value of the superionic cubic-superionic 
orthorhombic boundary slope at that triple point, and as- 
sumed it to be constant along the multi-phase boundary. 
(In fact, neglecting dT/dP variations of this type may 
introduce some errors in our predictions however, as we 
will show in Sec lIII Cl these turn out to be rather small.) 
By tracing the superionic fluoride-superionic PbCl2-type 
phase boundary, a second special triple point (-Pjc^l^c) 
is found at its intersection with the solid orthorhombic- 
superionic orthorhombic phase boundary (see Fig. [3]). 
Likewise, one-phase MD simulations were conducted at 
(P2c, 72c) to obtain the value of the corresponding dT/dP 
slope and assumed it to be constant along the superionic 
cubic-solid orthorhombic phase boundary. Proceeding 
as before, we identified the existence of a third special 
triple point (PscTsc) at the crossing of the superionic 
cubic-solid orthorhombic phase boundary with the solid 
fluoride-superionic fluoride phase boundary. Finally, the 
boundary separating the solid cubic and solid orthorhom- 



bic regions was drawn by joining the thermodynamic 
states (PscTsc) and {Pt,0), where Pt is the pressure at 
which the cubic — )■ orthorhombic transformation occurs 
at T = 0. 

In what follows we describe and yield the parametriza- 
tion of all the multi-phase boundaries cited above, ex- 
plaining the simulation procedures that we followed to 
obtain them. 



A. Superionicity 

Comprehensive one-phase {N, V, T) molecular dynam- 
ics simulations were carried out to compute the solid- 
superionic phase boundary of cubic and orthorhombic 
CaF2 as a function of pressure. Calculations comprised 
large simulation boxes of 6, 144 atoms and exceptionally 
long simulation times of ^ 200 ps. We systematically car- 
ried out simulations at temperature intervals of 100 K, 
up to 3000 K, at each volume considered. 

Following previous works 0, S HE H^j we identified 
superionicity via inspection of the time-dependent mean 
squared displacement function (MSD) obtained in one- 



6 



phase MD simulations. The MSD function is defined as 

{\ARnt)\)^{mt + to)-Mto)\') , (4) 

where Ri(t) is the position of atom i at time t, Iq is an 
arbitrary time origin, and (• • • ) denotes thermal aver- 
age. ((|Ai?|(t)|} was computed separately for each ionic 
species and thermal averages were performed over to and 
atoms -see Fig. HI-). In practice, F~ diffusion is signaled 
by the appearance of a non-zero MSD slope as we illus- 
trate in Fig. m It is worth noticing that anion diffusion 
is hardly discernable when scrutinizing only spatially av- 
eraged static quantities; for instance, one can not ap- 
preciate important differences between the CaF2 radial 
pair distribution functions obtained at temperatures just 
above and below the superionic transition (see Fig. [S]). 

For each phase, we determined the superionic transi- 
tion temperature Ts{P) at six different volumes. In both 
phases we found that the results could be perfectly fitted 
to straight lines Tg = as + bgP, where Os = 1379.69 K and 
bs — 34.23 K/GPa are the optimal parameters for the cu- 
bic phase and a, = 1497.70 K and bs = 50.17 K/GPa 
for the orthorhombic phase (see Fig. [3]). As it can 
be observed, the slope of both fluoride and PbCl2-type 
Ts{P) curves are positive so indicating that the supe- 
rionic phases are entropically stabilized over the corre- 
sponding crystals (i. e. the volume of the former systems 
is larger than those of the last -0 < AV- and consequently 
so is the enthalpy -0 < AH-; then AH = TAS must be 
accomplished over the coexistence line and < AS fol- 
lows). This result is reminiscent of customary melting 
however no depletion of the two-phase boundary induced 
by compression is observed in the present case. 

For ambient pressures, we obtain a superionic tran- 
sition temperature of Ts{0) — 1380 (10) K which is in 
very good agreement with ex per imental data T/^^'* — 
1420 (20) K found in Refs. [Mli. Also, we com- 
puted the diffusion coefficient of F~ anions, -D_, from 
a least squares fit to the MSD profiles (i.e. {\ARf{t)\) = 
A- + 6D-t, following the well known Einstein rela- 
tion). Our Z?_ results at temperatures 1400, 1500 and 
1600 K are 0.3, 0.9 and 2.7 10"^ cm^/s, which turn 
out to be consistent with experimental data 0.6, 1.6 and 
3.2 10-5 cmVs 0,i,|33,|33- 



B. Melting 



Following previous works [3J-|38[, we performed com- 
prehensive {N,V,E) two-phase coexistence MD simula- 
tions in order to determine the melting curve of cubic and 
orthorhombic CaF2 under pressure. Starting with a su- 
percell containing the perfect crystal structure (i.e. either 
fluoride or PbCl2-type), we thermalize it at a tempera- 
ture slightly below the expected melting temperature for 
about 10 ps. The system remains in a superionic state. 
The simulation is then halted and the positions of the 
atoms in one half of the supercell are held fixed while 
the other half is heated up to a very high temperature 
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FIG. 4: Mean squared displacement of fluor and calcium 
species represented as a function of time at temperatures be- 
low (Top) and above (Bottom) the corresponding transition 
point Ts — 1550 K. One-phase molecular dynamics simula- 
tions were performed for cubic CaF2 at pressures ~ 5.3 GPa. 



(typically five times the expected melting temperature) 
for about 60 ps, so that it melts completely. With the 
fixed atoms still fixed, the molten part is rethermalized 
to the expected melting temperature (for about 10 ps). 
Finally, the fixed atoms are released, thermal velocities 
are assigned, and the whole system is allowed to evolve 
freely at constant {N, V, E) for a long time (normally 
more than 100 ps), so that the solid and liquid come into 
equilibrium. The system is monitored by calculating the 
average number of particles in slices of the cell taken par- 
allel to the boundary between the solid and liquid. With 
this protocol, there is a certain amount of trial and er- 
ror to find the overall volume which yields the coexisting 
solid and liquid system. (An example of a successful co- 
existence run is shown in Fig. [6]) Our simulations were 
done on cells containing 12, 288 atoms with the long axis 
being perpendicular to the initial liquid-solid boundary. 

In both cubic and orthorhombic structures, we calcu- 
lated the melting transition temperature Tm at six dif- 
ferent volumes. We find that our results can be very well 
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FIG. 5: Radial pair distribution functions obtained at tem- 
peratures below (Top) and above (Bottom) the correspond- 
ing superionic transition point Ts — 1550 K of cubic CaF2 at 
P - 5.3 GPa. 

fitted to the so-called Simon equation [iOl 

Tr^{P) = a (^1 + , (5) 

where a = 2044.05 K, b = 1.2049 GPa and c = 0.1202 are 
the optimal values for the cubic phase, and a = 389.19 K, 
b = 0.0180 GPa and c = 0.3159 for the orthorhom- 
bic. Our predicted melting temperature at P = is 
2044 (100) K, which turns out to be slightly larger than 
the experimental value 1690 (20) K [111,53]. At this tem- 
perature the calculated bulk modulus of the superionic 
phase is ^ 57 GPa, which is significantly smaller than the 
obtained for the perfect cubic crystal (see Table II). Using 
the Simon equation, we estimate the corresponding zero- 
pressure melting slope to be dT,n/dP — 203.97 K GPa^^. 
We note that this quantity is considerably reduced un- 
der pressure, in constrast to the constant dTg/dP case 
reported in the previous section; at P = 5 GPa, for in- 
stance, we find dTm/dP '--^ 40 K GPa^^. Interestingly, 
the value of the calculated high-P melting slopes are 
larger than those measured in alkaline-earth (AE) flu- 
orides (e.g. LiF and NaF) 43], in spite of the structural 



FIG. 6: Determination of melting temperatures based on 
two- phase coexistence molecular dynamics simulations. Top: 
Number of particles histogram represented as a function of 
position along the direction parallel to the initial solid-liquid 
boundary. Bottom: Atomic positions projected over the a; — y 
plane of the simulation box. Calcium and fluor ions are rep- 
resented with red and green dots, respectively. 



similarities between both type of structures. AE-fluorides 
(rock-salt) and CaF2 (fluorite) have both an fee array of 
anions and differ only in the positions of the cations. A 
possible cause for the larger melting slope of GaF2 could 
be that AE-fluorides melt directly from an ordered crys- 
talline structure, while CaF2 melts from the superionic 
phase. Due to diffusion of F atoms, the degree of atomic 
disorder is larger in superionic CaF2 than in rock-salt 
AE-fiuorides. Therefore, it is smaller the corresponding 
entropy of fusion (and hence so AH). Consequently, ac- 
cording to the Clasius-Clayperon relation, the slope of 
the melting line should be larger in CaF2. Regarding the 
melting properties of CaF2 in the orthorhombic phase, 
this exhibits a very low zero-pressure melting point of 
^ 400 K (a temperature at which the corresponding crys- 
tal structure is metastable) and a overall steep dT^/dP 
slope (see Fig. It is worth noticing that the melting 
line of the orthorhombic phase intersects that of the cubic 
phase at the special triple point (7.18 GPa, 2580 K). At 
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the triple point an increase occurs in the melting slope. 
This phenomenon is most likely due to the volume change 
associated with the solid-solid transition Unfortu- 
nately, we do not know of any experimental data to com- 
pare with these results. 

Concerning previous estimations, we are just aware 
of works done by Zijiang [Tsl], Wang fl^, and Zeng et 
al. [l7j who employed similar pair-potential models than 
us here. The P = melting temperatures predicted by 
those authors range from 1650 to 2100 K, which are in 
reasonable good agreement with our results. Neverthe- 
less, we must note that the computational approaches 
used in those studies mainly consist of one-phase MD 
simulations and heuristic overheating arguments, both 
of which are well-known to produce very im prec ise T,„ 
results [Hil m Ei- In fact, Zeng et al. [l^ find a 
melting temperature of 990 — 1073 K for the Pnma phase 
of CaF2 at P = 10 GPa, a result that appreciably differs 
from ours Tm = 2867 K obtained under identical condi- 
tions. The main reasons for such discrepancies may prob- 
ably lie on the different methodologies employed and the 
complete omission of superionic effects in Zeng's work. 



C. Special triple points and other two-phase 
boundaries 

In previous sections, we have characterized four out of 
the seven two-phase boundaries shown in Fig. |31 As it 
has been explained at the beginning of this section, we 
determined the remainder of boundaries (i.e. solid-solid, 
solid-superionic and superionic-superionic) based on fun- 
damental thermodynamic considerations, one-phase MD 
simulations, and assuming linear pressure dependence in 
all of them. The results so obtained are summarized 
in Table IV. As one can observe, the slope of the su- 
perionic fluoride-superionic PbCl2 boundary is infinite 
because the calculated enthalpy difference between the 
two superionic phases at (7.18 GPa, 2580 K) is zero (see 
Eq. [3]). In the other two cases, we find that the slope of 
the superionic fiuoride-PbCl2 boundary is positive (i.e. 
AH < and AV < 0) and roughly a factor of two larger 
in absolute value than of the fluoride-PbCl2 boundary. 
In this last case, the resulting dT/dP slope is negative 
because the differences in enthalpy and volume between 
the two crystals are of opposite sign (i.e. AH > and 
AV < 0). In fact, a priori one would expect the entropy 
of the fluoride — >■ PbCl2 transformation to be positive, 
and hence so AH, because higher symmetry structures 
(i.e. cubic CaF2) in general imply lower entropy. 

Interestingly, we identify the presence of three spe- 
cial triple points, (PcTc), in the thermodynamic region 
6 < P < 8 GPa and 1500 < T < 2750 K (see Table IV 
and Fig. [3]). These special thermodynamic states are lo- 
cated at the intersections between three different phase 
boundaries, and the pressures and temperatures at which 
are predicted to occur in principle can be accessed in 
experiments. To this regard, information contained in 



Fig. [3] and Table IV must be considered as highly valu- 
able since identification of coexisting superionic and liq- 
uid phases turns out to be very challenging in practice. 
In fact, we just know of a couple of recent experimental 
works wherein coexistence between superionic and liquid 
phases has been suggested to happ en in water upon very 
extreme P — T conditions [13, l48j |. 

If one lifted the linear pressure dependence assump- 
tion from the solid-solid, solid-superionic and superionic- 
superionic phase boundaries, the location of the three 
special triple points (Pc,Pc) quoted in Table IV will 
probably change. In order to quantify the magnitude 
of those variations and to assess so the accuracy in our 
results, one could for instance perform calculations of 
the Gibbs-Duhem integration type and exactly determine 
the involved multi-phase boundaries |49l - [5l| . Gibbs- 
Duhem integration calculations and other equivalent ex- 
act schemes [H, [H, however are computationally 
very intensive so that we opted for a more straight- 
forward test. In particular, we computed the value of 
the dT/dP slope at states {P'^,T'^ found at halfway 
of the approximated linear multi-phase boundaries (see 
Fig. [3]), and checked whether these differed apprecia- 
bly or not from those reported in Table IV. For the 
superionic cubic-superionic orthorhombic phase bound- 
ary, we find that the enthalpy difference between the two 
phases at (7.18 GPa, 2220 K) is zero implying also an in- 
finite slope; assuming linear pressure dependence, there- 
fore, seems to be adequate in this case. By contrast, 
at point (6.98 GPa, 1728 K) belonging to the superionic 
fluoride-solid PbCl2 boundary we obtain a dT/dP value 
that is roughly twofold larger than the obtained at the 
corresponding special triple point (7.18 GPa, 1860 K) . 
However, the extend of this last two-phase boundary is 
so reduced that we may still assume that the resulting 
{Pc,Tc) inaccuracies are reasonably small. Certainly, at 
state (8.87 GPa, 800 K) of the cubic-orthorhombic phase 
boundary we find dT/dP = -300 (30) K/GPa which is 
in very good agreement with the constant assumed value 
of -315 (40) K/GPa reported in Table IV. These out- 
comes come to show that assuming linear pressure depen- 
dence in all two-phase boundaries involving superionic 
and crystal structures provides very consistent results, 
in spite of the small errors introduced in the superionic 
cubic-solid orthorhombic boundary. Consequently, our 
special triple point estimations can be safely considered 
as accurate. 



IV. THE ROLE OF REPULSIVE AND 
DISPERSION INTERACTIONS ON MELTING 

In contrast to ab initio electronic band structure meth- 
ods, empirical and semi-empirical force fields may suffer 
from versatility and transferability issues. This means 
that an interaction model which correctly describes a set 
of properties may fail at reproducing others and/or the 
same under different thermodynamic constraints. Actu- 
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1^ (K/GPa) Pc (GPa) Tc (K) AV (AVatom) AH (eV/atom) 



Fluoride - PhCh 


-315 (40) 


6.78 


(5) 


1600 


(10) 


-0.97 


(2) 


0.031 (5) 


Superionic Fluoride — PbCl2 


650 (100) 


7.18 


(5) 


1860 


(10) 


-1.01 


(2) 


-0.023 (5) 


Superionic Fluoride — Superionic PbCl2 


oo 


7.18 


(5) 


2580 


(10) 


-0.36 


(2) 


0.00 (5) 



TABLE IV: Thermodynamic data describing the unravelled solid-solid, solid-superionic, and superionic-superionic phase bound- 
aries. The slope of the P — T boundaries and corresponding changes in volume {V) and enthalpy (H) at the three special triple 
points {Pc,Tc) are reported. 



ally, a considerable number of CaF2 pairwise potentials 
are found in the literature each having been designed for 
a distinct purpose [7,, 9, .20, 52, 53j]. Trying to derive 
the phase diagram for all of them would indeed be a te- 
dious and extremely boring task. Fortunately, once the 
phase stability properties of a given interatomic potential 
are known it is possible to deduce those for other simi- 
lar interaction models in a computationally efficient and 
physically insightful way. We refer here to the original 
free-energy perturbative approach developed by Gillan 
and collaborators and which has been successfully ap- 
plied to the study of transition metals under extreme 
P — T conditions [s^ - IsgI [s^ [39|. In this section, we use 
Gillan's ideas to compute the shift in melting tempera- 
ture caused by mild variations of the potential parame- 
ters employed through this work (i.e. under the general 
transformation {Xij where {Xij} correspond 

to the set of parameters reported in Table I). The mo- 
tivation for this analysis is not only to gain access to 
the phase diagram features of other similar BMH poten- 
tials but to understand also the general role of short- and 
long-ranged interactions in melting. A brief description 
of Gillan's free-energy perturbative approach is provided 
next. 

For a given P and T, the difference Gjf = Gj, - 
between the Gibbs free energies of the new liquid and 
solid (i.e. obtained with the new set of potential param- 
eters {Xlj}) deviates from the corresponding difference 
Gq = Gq — Gq of the initial reference liquid and solid (i.e. 
obtained with the reference set of potential parameters 
{Xij} shown in Table I), and we write: 

Gjf (F, T) = G[," (P, T) -I- AG'" (P, T) . (6) 

The shift AG'''(P, T) caused by changing the total-energy 
function from Uq to Un, induces a shift in the correspond- 
ing mel ting temperature Tm{P). To first order, the latter 
shift is p39| : 

AT -^^^liZ^ (7) 

^-l-m— qls ' ) 

where Sq^ is the difference between the entropies of the 
liquid and solid (i.e. the entropy of fusion) of the initial 



reference system, and AG is evaluated at the melting 
temperature of the initial reference system. The shift 
AG'" is the difference of shifts of Gibbs free energies of 
the liquid and solid caused by the shift A[/ = [/„ — J/q- 
Under constant volume and temperature, the shift of 
Helmholtz free energy AF arising from AU is given by 
the well-known expansion: 

AF = (AC/)o - i/3(5AC/2)o + • • • , (8) 

where ^ = l/fceP, 5AU = AU- {AU)o, and the averages 
are taken in the initial reference ensemble. From AF, we 
obtain the shift of Gibbs free energy at constant pressure 
as: 

AG = AP - ^Fkt(AP)2 , (9) 

where kt stands for the isothermal compressibility and 
AP is the change of pressure caused by the replacement 
Uo Un at constant V and T. 

Our ATm calculations were done over series of 
AXij/Xij points {AXij = X[j — X^, where Xij refer to 
the reference potential parameters reported in Table I) 
generated in the range —25 % < AXij/Xij < 25 % and 
taken at 1 % intervals (see Fig.[7]). The averages involved 
in these calculations were computed over 500 liquid and 
solid configurations generated in long one-phase MD sim- 
ulations performed with the reference BMH potential. 
These MD simulations were carried out at the arbitrarily 
selected state (1.0 GPa, 2180 K) lying over the superionic 
fluoride-liquid phase boundary. 

Fig. [7] shows our AT™ results expressed as a function 
of AAij / Aij and Adj/Cij variations (i.e. we have ne- 
glected pij fluctuations). As one may appreciate, signifi- 
cant variations of the short-ranged and long-ranged parts 
of the interatomic potential in general have only a mod- 
erate effect on Tm- For instance, Aij and Cij relative 
variations of 25 % provoke at most melting tempera- 
ture shifts of ~ 100 K. This result is surprising but at 
the same time reassuring in the sense that it adds ro- 
bustness to our Tm conclusions drawn in the previous 
section. Now, let us focus on the ATm shifts caused by 
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FIG. 7: Shift of the melting transition temperature ATm 
caused by the variation of potential parameters at Pm ~ 
1.0 GPa and TW, = 2180 K. Top: Short-range repulsive Ach-f 
and ^F-F parameters are varied and the rest of parameters 
kept fixed. Bottom: Long-range attractive Cca-F and Cf-f 
parameters are varied and the rest of parameters kept fixed. 



individual Ap-p, ^caF, C'ff, and CcaF variations. We 
note that due to the perturbative character of our ap- 
proach the total melting temperature shift provoked by 
a general {Xij -^Ij} transformation is equal to the 
sum of {ATm} shifts caused by the individual Xij fluc- 
tuations. Concerning the short-ranged part of the inter- 
action model (see top of Fig. [7]) , we find that reduction 
of the ^FF parameter causes a positive melting temper- 
ature shift; by constrast, increase of the same parameter 
leads to AT^ < 0. This outcome comes to show that 
strengthening (weakening) of the repulsive F-F interac- 
tions tends to further stabilize (destabilize) the liquid 
phase. Interestingly, we observe the opposite trend for 
AcaF ■ AAcaF < fluctuatious lead to AT^ < (i.e. the 
liquid phase is energetically favored) whereas AAcaF > 
lead to AT„i > (i.e. the liquid phase is energetically 
disfavored). Upon a same AAij/Aij change, we observe 
that the melting temperature shift obtained in the Ap-p 
case is smallest (in absolute value) (e.g. AT^ = — 10 
and 51 K for AAffMff and AAcaFMcaF = 25 %, re- 



spectively) thus repulsive Ca-F interactions play a more 
dominant role in melting. Regarding the long-ranged 
attractive part of the interaction model (see bottom of 
Fig. EI), the situation is the opposite than just explained. 
In particular, ACff < (ACff > 0) changes lead to 
AT™ < (AT™ > 0) and ACcaF < (ACcaF > 0) to 
ATm > (AT„ < 0). Moreover, upon a same AC^j / C^j 
variation the melting temperature shift obtained in the 
Cff case is largest (in absolute value) (e.g. AT,„ — 83 
and -29 K for ACff/Cff and ACcaF/CcaF = 25 %, re- 
spectively) thus attractive F-F interactions play a more 
important role in melting. As a summary of these re- 
sults, we can state that short-ranged repulsive Ca-F and 
long-ranged attractive F-F contributions to melting are 
most notorious and that reduction (increase) of the in- 
volved paremeters Aq^f and Cff leads to further sta- 
bilization (destabilization) of the liquid over the cubic 
superionic phase. 

Finally, we can use our results shown in Fig. [7] to 
predict, at least at a qualitative level, the melting fea- 
tures corresponding to other CaF2 BMH potentials. In 
particular, we examine two parametrizations indepen- 
dently proposed by Gillan Q and Boulfelfel [52]. For 
the sake of simplicity, we consider here only ATm con- 
tributions stemming from the dominant ^caF and Cff 
parameters. For Gillan's parametrization, we obtain 
A^CaFMcaF 60 % and ACff/Cff 550 %. Actu- 
ally, these values turn out to be exceedingly large so as 
to being treated within our perturbative approach how- 
ever we can make a qualitative statement based on the 
size and sign of those deviations. In particular, the Cff 
difference is eminently the largest and positive so that 
the resulting temperature correction is very likely to be 
positive and large (i.e. Tm ^ 2180 K). In the case of 
Boulfelfel's potential, we obtain AAcaFMcaF 72 % 
and ACff/ Cff ~ —100 %. These values are of opposite 
sign to Gillan's and still too large so as to being ana- 
lyzed with our method. Nevertheless, at the qualitative 
level, we can state that the resulting AT„j difference is 
very likely to be negative and small since the Cff devia- 
tion is negative and only slightly superior than j4caF (i-e. 
T,„ < 2180 K). 

As a concluding remark to this section we want to men- 
tion that the computational strategy just presented can 
also be applied to the analysis of multi-phase boundaries 
others than melting, and in general to the modelling of 
atomic interactions for derivation of phase diagrams. 



V. CONCLUSIONS 

To summarize, we have studied the phase diagram 
of CaF2 under pressure using classical atomistic simu- 
lations and a simple pairwise interatomic potential of 
the Born-Mayer-Huggings form. Our results show that 
a rich variety of crystal, superionic and liquid phases co- 
exist within the thermodynamic region < P < 20 GPa 
and < T < 4000 K. In particular, we find seven dif- 
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ferent two-phase boundaries for all of which we provide 
an accurate parametrization. Interestingly, three special 
triple points are predicted to exist within the narrow 
and experimentally accessible thermodynamic range of 
6 < P < 8 GPa and 1500 < T < 2750 K. Indeed, we 
believe that these stimulating findings should encourage 
new experimental searches in CaF2 under elevated P — T 
conditions. Also, we have analyzed the role of short- 
ranged repulsive and long-ranged attractive atomic in- 
teractions in the prediction of melting points, with the 
finding that repulsive Ca-F and attractive F-F contribu- 
tions are most notorious. In order to get rid of possible 
versatility and transferability force-field issues, it would 



be very much desirable to conduct ab initio simulation 
studies similar to the one presented here. Work in this 
direction is already in progress within our group. 
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